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Abetraot. 

This  paper  describes  a  solution  to  the  flexural  analysis 
of  a  beam  by  digital  oomputer.  Subroutines  were  written 
which  calculate  shear,  moment,  slope  and  deflection  at  any 
specified  point  on  a  beam  if  the  loading  is  known.  Shear 
defleotion  may  be  included  for  many  support  conditions. 
Neither  elastic  supports  nor  column  effeots  are  considered. 
The  subroutines  solve  problems  of  variable  beam  oross-seotion 
and  loading  by  utilizing  an  extension  of  "McCaulay's  Method", 
a  generalized  step  funotion. 

The  main  program  provides  input-output  facilities  and 
also  solves  for  indeterminate  reaotions  on  the  beam. 

Problems  and  their  solutions,  showing  the  versatility 
of  this  program,  are  also  inoluded. 
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1.  Introduction. 

The  routine  solution  to  engineering  problems  by  digital 
computer  appears  to  be  limited  only  by  the  scope  of  the 
available  programs.  A  flexible  program  permitting  the  solu¬ 
tion  of  routine  problems  in  flexural  analysis  of  beams  would 
seemingly  fill  a  need  in  the  field  of  digital  computer  ap¬ 
plications  to  Mechanics.  The  present  work  is  devoted  to  suoh 
an  idea. 

Five  subroutines,  "LOAD",  "SHEAR",  "MOMENT",  "SLOPE", 
and  "DEFLECT",  were  written  which  will  calculate  the  load, 
shear,  moment,  slope,  and  deflection,  respectively,  at  any 
indicated  point  on  a  beam  if  the  indeterminate  quantities 
are  calculated  elsewhere.  The  use  of  generalized  step  functions 
in  this  project  allows  for  piecewise  linear  variation  of 
distributed  loading,  bending  compliance,  and  shear  com¬ 

pliance, 

The  main  program,  "BEAM3",  first  generates  and  then 
solves  a  set  of  simultaneous  linear  equations,  calculating 
the  indeterminate  quantities  of  the  beam.  The  above  mention¬ 
ed  subroutines  are  then  utilized  to  calculate  the  remaining 
flexural  quantities.  This  program  also  provides  input-output 
facility. 
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2.  Definition  of  the  problem. 

The  problem  of  structural  beam  analysis  is  one  in 
which  the  external  loading  applied  to  the  beam  is  well  de¬ 
fined,  as  is  the  geometric  configuration  and  the  nature  of 
the  beam  itself.  Where  information  about  such  a  beam  is  re¬ 
quired,  this  information  will  most  likely  fall  into  the 
following  categories: 

1.  Internal  shear  forces,  V,  and  moments,  M. 

2.  Slope,©. 

3.  Deflection,  y. 

4.  Redundant  reactions,  either  forces,  Ry,  or  moments, 

RMZ. 

"Flexural  analysis"  is  used  here  as  meaning  the  process  by 
which  this  information  is  generated  from  the  known  data. 

In  attacking  this  problem  of  flexural  analysis,  the 
intention  was  to  limit  the  scope  as  little  as  possible. 
However,  several  customary  simplifications  and  limitations 
were  made.  First,  elastic  supports  and  beam  column  effects 
were  not  considered.  Limitations  concerning  the  type  and 
number  of  external  loads  and  the  geometry  of  the  beam  will 
be  discussed  later.  Also  the  basic  relations  which  follow, 
inherently  restrict  the  deflection  of  the  beam  to  small  values 
and  the  internal  stresses  to  values  below  the  elastic  limit 
of  the  structural  material.  Shear  deflection,  which  is 
generally  neglected  in  a  problem  of  this  type,  has  been 
Included  in  some  problems, 
a.  Bending  deflection  (y^) . 

The  Bernoulli-Euler-Navier  equation  governing  the  elastic 
curve  of  a  beam  due  to  bending,  consistent  with  the  sign 
convention  shown  in  Appendix  B  is: 
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«  curvature  *  (d2y/dx2)  [l+(dy/dx)2]  ^ 

2  2 

and  we  use  the  approximation,  curvature  d  y/dx  . 

The  following  equations,  derived  from  the  above  equation 
and  shown  in  Timoshenko  /5/,^  with  differences  due  to  the  chosen 
sign  convention,  describe  the  flexure  of  a  beam  (See  Appendix  B 
for  notation). 


W  *  q(x)  , 
V  *  f\ Wdx  , 
M  *  _/Vdx  , 


yi  -  Jr91  dx  , 
b.  Shear  deflection  (y2). 

The  equation  governing  the  shear  deflection  of  a  beam, 

also  given  in  Timoshenko  /5/»  compatible  with  the  stated 

sign  convention  becomes: 

dy~  KV 

dx  AG 

It  is  evident  that: 


o.  Total  deflection  (y). 

The  equations  used  in  this  thesis  are  the  result  of  the 
sum  of  the  above  shear  and  bending  deflections.  They  are: 

W  =  q(x)  , 


Numbers  so  indicated  refer  to  the  bibliography. 
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V 


M 

e 

Y 


+  EFj* 

+ 

KV 

AG 

♦  c2** 


♦ 


Note  that  step-discontinuities  may  appear  in  all  of  the 
above  equations,  except  in  the  equation  for  deflection.  These 
discontinuities  limit  the  methods  that  can  be  used  in  handling 
the  equations. 

Three  methods  that  could  be  used  to  solve  the  equations 
are  a  pure  analytical  approach,  a  step-by-step  numerical 
integration  scheme,  and  the  "curly  bracket"  method.  Of 
these  methods,  the  analytical  approach  is  generally  restrict¬ 
ed  to  a  specific  problem  and  is  unsuited  for  adaptation  to 
digital  computer  programming.  The  numerical  integration  scheme 
can  be  used  in  the  computer  solution  of  flexural  analysis 
and  may  have  been  programmed;  however,  our  search  of  literature 
has  not  uncovered  any  such  program.  The  "curly  bracket" 
method  was  chosen  for  use  in  this  thesis  because  it  was  felt 
that  it  would  result  in  a  more  exact  solution  requiring  less 

*In  these  equations,  2^  and£M  represent  the  sum  of 
the  point  forces  and  moments,  respectively,  applied  to  the 
beam. 

**C^  and  C2  are  constants  of  Integration. 
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3.  "Curly  bracket"  method 

The  "curly  bracket"  method  is  based  upon  the  properties 
of  a  convenient  notation  defined  as  follows: 

{x-a} n  =  0  i  f  x  <a 

=  ( x- a ) n  if  x> a 

where  n  is  an  integer^O. 

Appendix  A  contains  further  properties  of  expressions  contain¬ 
ing  curly  brackets.  The  notation  is  generally  credited  to 
McCaulay  /l/,  but  similar,  though  less  useful  notations  can 
be  traced  to  A.Fdppl  and  St.  Venant.  Also,  the  ideas  in¬ 
volved  are  related  to  the  use  of  the  Laplace  Transformation 
with  the  result  that  this  method  or  approximations  to  it  are 
reinvented  frequently  /&/  /9/. 

Using  this  method  allows  one  equation,  valid  over  the 
entire  length  of  the  beam,  to  be  expressed  for  each  of  the 
basic  equations  of  flexure.  The  particular  value  of  this 
method  is  that  it  permits  the  integration  of  these  equations 
to  be  carried  out  easily  and  has  the  virtue  of  adjusting  the 
constants  of  integration  most  conveniently. 

In  addition,  for  use  with  a  digital  computer,  the  "switch¬ 
ing"  property  of  the  curly  bracket  can  be  easily  handled  in 
FORTRAN  by  use  of  an  "IF"  statement. 

A  simple  example  should  convince  most  readers  of  the 
convenience  of  this  method.  Examine  the  following  problem, 
which  consists  of  a  uniform  beam,  simply  supported  at  each 
end  with  a  concentrated  load  at  the  center  of  the  span. 
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stant  variation  of  bending  compliance,  In  this  thesis 

the  method  has  been  further  expanded  to  include  beams  with 
piecewise  linear  variation  of  both  bending  compliance  and 

,  K 

shear  compliance,  . 

AG 
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4.  Implementation. 

Although  the  formulas  shown  in  Appendix  A  could  be  ex¬ 
panded  further,  a  piecewise  linear  variation  of  distributed 
loads,  bending  compliance,  and  shear  compliance  was  consider¬ 
ed  capable  of  providing  sufficient  accuracy  and  flexibility 
for  the  resulting  program. 

Fig.  2(a)  Illustrates  a  distributed  load  on  a  beam 
whose  origin  is  at  x=0.  Fig.  2(b)  shows  an  element  of  the 
distributed  load. 


Fig.  2(b) 


The  analytical  expression  for  this  element  is: 


Q(x)  -  QaJx-A}0  +  ^^(x-A}1  -  ^{x-B]0  -  {x-B> 1 


The  same  relationship  Is  used  for  an  element,  f(x),  of  bend¬ 
ing  compliance  and  an  element,  g(x),  of  shear  compliance  upon 
substituting  f  for  Q,  and  g  for  Q,  respectively. 

The  basic  equations  of  flexure  can  now  be  expressed  In 
the  same  manner  as  they  were  utilized  In  programming.  They 
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are; 


NOL 

W(x)  -  2  Q,(x) 

1=1  1 


(1) 


£NOF  .  NOR 

W(x)dx  ♦  2 


(2) 


where  F  are  upward,  known  forces  at  x^x*.  and  R  are  up- 
yl  1  yi 

ward,  unknown  reactions  at  x^x 


NOM 


NORM 


M(x)  =  \  V(x)dx  +  2  fx~x  )°  +  RM  £x-x  V 

Jo  i=l  zi^  V  i-1  zi*  rmii 

(3) 

where  M  are  clockwise,  known  moments  at  x-x  and  RM,  are 
zi  mi  zi 

clockwise,  unknown  reactive  moments  at  x«x 


rnu 


G 


4y  = 

dx 


NOI 


M(x)  X!  f*(x)dx  - 
i=l 


NOK 

v(x)  X  e«(*} 

i=l  1  . 


♦  c- 


(4) 


where  is  the  constant  of  integration. 

Appendix  A  indicates  how  the  product  and  the  integration  of 
the  product  of  the  curly  bracket  symbols  are  formed. 


^xrxr  _  NOI 

M( 


inuj.  -1  ,-x-  NOK  -I 

1)2  fi(  %  )J  dx  v(x)  2g1(x)dx  +  +  C2 

i-1  -  1-1  (5) 

where  C 2  is  the  constant  of  integration. 

Equations  (1)  through  (5)  when  properly  coded  in  FORTRAN 


The  numbers,  NOL,  NOF,  NOR,  NORM,  NOI  and  NOK,  appear¬ 
ing  above  the  summation  signs  are  program  input  parameters 
indicating  the  number  of  terms  for  each  problem. 
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form  the  core  of  the  five  basic  subroutines?  LOAD,  SHEAR 
MOMENT,  SLOPE  and  DEFLECT.  It  must  be  noted  before  discuss¬ 
ing  the  subroutines  that  equations  (2)  through  (5)  Involve 
the  unknown  values  C2,  Ry  ,  and  RMz  ,  where  there  may  be 
as  many  as  ten  values  each  of  Ry  and  RMZ.  The  calculation  of 
these  values  is  done  in  the  main  program  BEAM3  and  will  be 
discussed  in  section  6. 
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5.  The  basic  subroutines 


Subroutines  LOAD,  SHEAR  and  MOMENT, 

The  above  named  subroutines,  as  their  names  imply,  are 
employed  to  calculate  the  total  distributed  load,  shear  and 
moment,  respectively,  at  any  desired  point  on  the  beam.  The 
expanded  forms  of  equations  (1),  (2)  and  (3)  were  used  in 
writing  these  subroutines.  Utilizing  the  distance,  x^,  as 
the  entering  argument,  the  subroutines  select  and  sum  the 
appropriate  terms  in  the  corresponding  equation  by  inspect¬ 
ing  all  loading  terms  and  discarding  those  in  which  the  value 
within  the  curly  bracket  is  negative.  Actually  the  load, 
shear  and  moment  are  calculated  a  small  increment,  £  -  10  , 

to  the  right  of  x^  in  order  to  obtain  the  right  side  value 
of  any  step-discontinuities  occurlng  at  the  point,  x ^  (See 
flow  charts  on  pages  13,  14  and  15). 

Subroutine  SLOPE. 

This  subroutine  is  equation  (4),  excluding  the  constant 
of  integration,  C^,  in  FORTRAN  language  (See  flow  chart  on 
page  17).  Note  that  the  equation  can  be  expressed  as  a 
summation  of  terms,  which  are  constants  times  the  multiplica¬ 
tion  of  two  step- functions  of  the  form  C  £x-aj  n£x-b} m,  or  the 
summation  of  the  integral  of  the  same  type  of  terms,  where  n 
equals  either  0  or  1. 

Function  subroutines,  FTON  and  FT1N  (see  Appendix  A), 
will  calculate  the  value  of  £x-a}  °  [x-bjm  and  [x-a^  ^{x-b}  n, 
respectively.  Function  subroutines,  ENTON  and  ENT1N  (See 
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Flow  Chart  of  LOAD 


Note:  Numbers  preceded  by  a  number  sign,  refer  to 

the  statement  number  in  the  listing  of  LOAD. 
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Flow  Chart  of  SHEAR 


Note:  Numbers  oreceded  by  a  number  sign,  refer  to 

the  statement  number  In  the  listing  of  SHEAR. 
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Flow  Chart  of  MOMENT 


RETURN 

Note:  Numbers  preceded  by  a  number  atgr,  #,  refer 
statement  number  In  the  listing  of  MOMENT, 
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Appendix  A),  will  calculate  the  first  integral  of  these 
functions.  Equation  (4)  can  now  be  expressed  as  the  summation 
of  constants  multiplied  by  the  appropriate  function  subroutines. 

SLOPE  calculates  all  combinations  of  M(x)  and  f(x)  by 
stepping  through  all  the  f(x)  terms  and  varying  all  the  terms 
in  the  moment  equation  at  each  step.  It  then  steps  through 
all  the  g(x)  terms  varying  all  the  terms  in  the  equation  for 
shear  at  each  step.  The  sum  of  these  terms  is  the  value  of 
the  slope,  minus  the  integration  constant,  0^,  at  the  entered 
distance,  x. 

Subroutine  DEFLECT. 

At  this  point  the  author  wishes  to  point  out  the  similar¬ 
ity  of  the  equations  for  slope  and  deflection.  Since  the 
integration  constants,  C-^  and  C2,  are  again  excluded  from  the 
subroutine,  and  the  function  subroutines,  DITON  and  DITIN 
(See  Appendix  A),  calculate  the  seoond  integral  of  the  mul¬ 
tiple  step- functions,  subroutine  DEFLECT  is  identical  to 
SLOPE  except  for  substituting  ENTON  for  FTON,  ENT1N  for  FT1N, 
DITON  for  ENTON  and  DITIN  for  ENT1N.  The  result  is  a  sub¬ 
routine  that  will  calculate  the  deflection,  minus  the  con¬ 
stants  of  integration,  at  any  required  point  on  the  beam. 
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6.  Program  BEAM3. 

The  main  program,  BEAM3,  can  be  divided  into  four  dis¬ 
tinct  steps.  These  steps  are:  reading  the  input  data,  cal¬ 
culation  of  indeterminate  quantities,  calculation  of  the 
flexural  quantities,  and  the  output  of  the  solution. 

The  program  first  reads  the  input  data  which  includes 
the  length  of  the  beam,  all  external  loading,  the  location  of 
all  reactions,  the  flexural  properties  of  the  beam,  and  a  group 
of  parameters  giving  the  number  of  each  type  of  load  or  re¬ 
action.  It  also  provides  for  reading  a  predetermined  value 
of  deflection,  PDEFT,  at  each  reaction  and  a  predetermined 
value  of  slope,  PSLP,  at  each  moment  reaction. 

The  program  then  proceeds  to  solve  for  the  indeterminate 
quantities  of  the  beam  by  utilizing  the  two  equations  of 
static  equilibrium  plus  equations  obtained  from  the  known 
boundary  condition  at  each  reaction  or  reactive  moment.  The 
first  two  equa;ions  are  simply  the  summation  of  the  forces 
on  the  beam  equals  zero  and  the  summation  of  the  moments 
about  the  right  end  of  the  beam  equals  zero.  The  remainder 
of  the  equations  are  obtained  from  the  fact  that  the  slope 
at  each  reactive  moment  or  the  deflection  at  each  reaction 
is  known.  Using  this  fact  in  the  appropriate  one  of  either 
equation  (4)  or  (5)  results  in  a  set  of  linear  equations  for 
the  indeterminate  quantities. 

The  equations  thus  obtained  are  generated  in  a  matrix 
form  that  is  compatible  with  subroutine  GAUSS2  (See  Appendix  E) . 
This  subroutine  is  then  called  to  solve  the  equations  simulta- 
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neously.  If  the  matrix  is  singular,  GAUSS2  exits  to  the  call¬ 
ing  program  which  in  turn  prints  "MATRIX  SINGULAR"  and  then 
stops.  If  not,  the  subroutine  solves  the  equations  and  re¬ 
turns  the  calculated  values  to  the  main  program.  The  most 
serious  limitations  of  this  program  are  imposed  by  the  gen¬ 
eration  or  solution  of  the  simultaneous  equations,  therefore 
it  is  appropriate  to  include  and  explain  several  of  them  at 
this  point. 

First,  at  least  one  reaction,  even  if  it  is  later  forced 
to  equal  zero,  must  be  included  with  the  input,  or  a  singular 
matrix  will  be  generated. 

Including  shear  deflection  in  the  solution  must  also  be 
greatly  limited  at  this  point.  Since  in  general  the  shear 
includes  step  discontinuities,  they  will  also  appear  in  the 
slope  if  shear  deflection  is  included.  The  program  generates 
an  equation  for  slope  at  each  moment  reaction.  If  by  chance, 
either  a  reaction  or  point  force  occurs  at  the  same  point  as 
a  moment  reaction,  the  slope  at  this  point  is  double  valued. 
The  program  is  unable  to  distinguish  which  of  the  two  values 
is  required  except  when  it  occurs  at  the  extreme  ends  of  the 
beam.  For  example  a  beam  "built  in"  at  both  ends  is  accept¬ 
able.  If  the  beam  is  "built  in"  at  only  one  end,  it  must  be 
solved  in  a  manner  such  that  the  "built  in"  is  at  the  origin 
of  the  beam. 

The  program  then  divides  the  beam  into  a  number  of 
increments,  NOP-1,*  and  calls  the  subroutines,  LOAD,  SHEAR, 

*NOP  is  another  parameter  read  in  at  the  start  of  BEAM3, 
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MOMENT,  SLOPE  and  DEFLECT  to  calculate  the  flexural  quantities 
at  each  of  these  increments. 

The  output  consists  of  printing  the  distance  from  the 
left  end  of  the  beam  to  the  point  in  question  and  the  values 
of  the  load,  shear,  moment,  slope  and  deflection  at  that 
point.  The  values  of  all  of  the  indeterminate  quantities  are 
also  printed  out. 
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Fig*  7  Simplified  Flow  Chart  of  BEAM3 
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7.  Testing  of  Program. 

The  testing  of  the  program  and  the  subroutines  could, 
in  theory,  continue  indefinitely  because  of  the  unlimited 
number  of  beam  configurations.  The  author  has  attempted  to 
show  (See  Appendix  D)  solutions  to  test  problems  of  simple 
configurations  and  constant  bending  and  shear  compliance  to 
compare  with  classical  solutions.  Then  a  statically  deter¬ 
minate  problem  with  a  varying  bending  compliance  was  solved 
and  compared  to  a  solution  obtained  from  a  numerical  in¬ 
tegration  scheme.  The  most  rigid  test  for  the  program  was 
in  constructing  a  problem  for  a  beam  of  unlikely  configuration 
(See  Problems  6  and  7  of  Appendix  D),  and  solving  the  flexural 
analysis  of  this  beam  from  both  ends.  Note  that  the  only 
significant  difference  between  the  two  solutions  is  the  opposite 
signs  for  shear  and  slope,  which  would  of  course,  be  expected. 
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8.  Conclusions  and  Recommendations. 

The  thesis  shows  a  method  of  including  piecewise  linear 
variation  of  both  bending  compliance  and  shear  compliance  in 
the  flexural  analysis  of  beams  by  extending  "McCaulay's  Method". 
Adapting  this  idea  to  the  digital  computer  resulted  in  a 
program  capable  of  analyzing  a  very  general  type  of  beam. 
However,  had  time  permitted,  several  changes  would  have  been 
made  to  the  program.  First,  internal  monitoring  of  the  input 
data  would  be  an  asset  to  any  user.  Also,  using  the  same 
general  method  of  attack,  this  program  would  have  been  ex¬ 
panded  to  include  bending  compliance  and  shear  compliance  of 
a  piecewise  parabolic  nature. 

Other  additions  to  the  program,  much  larger  in  scope, 
would  be  including  elastic  supports  and/or  beam  column  effect. 

It  would  be  interesting  to  compare  the  results  of  this 
program  to  those  obtained  from  a  program  utilizing  the  numer¬ 
ical  integration  method.  Also,  a  more  exhaustive  study  of  the 
program  could  be  made  to  determine  accuracy  of,  and  time 
requirement  for,  solutions. 
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APPENDIX  A 


"Curly  bracket"  method,  definitions  and  formulas 
We  define: 

1.  [x-a]n  =  (x-a)n  ,  x >  a 
=  0  ,  x  <  a 

where  n  is  an  integer^  0 

From  this  definition,  it  is  possible  to  establish  the 
following  formulas.  Although  it  would  be  possible  to  express 
these  results  in  terms  of  curly  brackets,  in  most  cases  the 
resulting  expressions  would  be  rather  complicated.  The  forms 
exhibited  below  have  the  advantage  of  being  most  readily 
conformable  to  the  decision  commands  available  in  FORTRAN. 


(x-a) 

n+1 


n+1 


0 


x>  a 

Otherwise 


3.  [x-a]0;x-b3n 

4.  ^x-aj-^x-b}11 

5.  lX[x-a30[x-b}ndx 


( x-b ) n 
0 


x>  a  and  x>b 
Otherwise 


(x-b)n+1  +  (b-a)(x-b)n  ,  x  >  a 

and  x  >b 


(x-b) 

n+1 


A*=*l 

n+1 


n+1 


Otherwise 
n+1 


(a-b) 

n+1 


n+1 


x  >a  >b 


x>b>a 


0 


,  Otherwise 
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6.  f  £x-a}  1  [x-bj  ndx 
Jo 


(x-b) 

n+2 


n+2 


,  ^\nf2 

(a-b) 

(n+1) (n+2) 

=  (x-b)  + 

n+2 


+  (b-a)  (x-b) 

n+1 

,  x>a>b 
(b-a) (x-b) n+1 


n+1 


n+1 


x>b>a 


Otherwi se 


7-  f  f 

(x-b)"*2 

(n+1 


n+1) (n+2) 

,  ,  .n+2 

(a-b) _ 

(n+1) (n+2) 

,  v  n+2 
(x-b) 

(n+1) (n+2) 
0 


(x-a) (a-b) 

(n+1) 


n+1 

x  >a  >b 

x  >b  ;>a 
Otherwise 


»•  a  ^-*3i^-b}n^ta 


(x-b) 


n+3 


(n+2) (n+3) 


{  —a,  (a-b) 
+  (  } )  n+2) 

,  •u\rH’3 

(x-b) 

(n+2) (n+3) 


n+2 


n+2 


(b-a) (x-b) 
(n+1) (n+2) 

2<a-b)n*3 


(n+1) (n+2) (n+3)  , 

n+2 


(b-a) (x-b) 
(n+1) (n+2) 


x>a  >b 
x>  b  >a 


,  Otherwise 


Function  Subroutines  FTON,  FT1N,  ENTON,  ENT1N,  DITON,  DITIN 
The  above  named  function  subroutines  were  written  and 


utilized  to  evaluate  each  of  the  preceeding  formulas  numbered 
3  through  8.  The  values  x,  a,  b  and  n  were  used  as  entering 
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arguments  to  these  functions.  The  function  then  became  the 
variable  as  follows: 


FTON 

=  $x-  a]°ix-b}n 

FT  IN 

*  [x-a}  1  £x-bj  n 

ENTON 

=  C  [x-e]  °  [x-b]  ndx 

SO 

ENT1N 

X 

G 

S' 

1 

H 

or 

I 

& 

II 

DITON 

- ,  Go  -aJ°  ^-bina?  ^ 

DITIN 

=  ££2H1&-6}“4*te 

As  examples,  flow  charts  for  FTON  and  ENTON  are  shown  in 
Figs.  8  and  9. 
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Fig.  8 


Flow  Chart  of  FTON 
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Fig.  9  Flow  Chart  of  ENTON 
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APPENDIX  B 


Sign  Convention  and  Notation 

Sign  Convention 

The  following  sign  convention  has  been  employed: 

y 

t 

- i-  x  y  *  deflection 


O  =  slope  = 

dx 


»(r  i) 
v  tczzu 


M  »  bending  moment 


V  =  shear  force 


W 


TTTTT 


¥  =  distributed  loads 


Note:  1.  Point  forces,  Fy,  and  point  reactions,  R  ,  are 

positive  in  the  positive  y  direction.  y 

2,  Point  moments,  Ma,  and  point  reactive  moments, 
RM2,  are  positive  in  a  clockwise  sense. 


Notation: 

THESIS 

W 

V 

M 

e 

y 


PROGRAM 

C 

W 

V 

AM 

Y1 

Y 


Length  of  beam 

Total  distributed  load  at  a  point 

Shear 

Moment 

♦ 

Slope 

Deflection 
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Fy 

FY 

Point  load  or  force 

Ry 

RY 

Point  reaction 

Mz 

AMZ 

Point  moment 

RMZ 

RAMZ 

Point  reactive  moment 

Q(x) 

Element  of  distributed  load 

Qa 

AQ 

Distributed  load  (left  end) 

% 

QjB 

Distributed  load  (right  end) 

f(x) 

EA 

Element  of  bending  compliance 

(It) 

Bending  compliance 

(left  end  of  element) 

EB 

Bending  compliance 
fright  en<?  of  element) 

g(x) 

Element  of  shear  compliance 

(K_) 

VAG; 

GA 

Shear  compliance 
(left  end  of  element) 

GB 

Shear  compliance 
(right  end  of  element) 

Xf 

XFO 

Distance  to  Fy 

xr 

XR 

Distance  to  Ry 

Xm 

XM 

Distance  to  Mz 

xm 

RXM 

Distance  to  RMZ 

Ag 

AG 

Distance  to  Ga 

Bg 

BG 

Distance  to  G^ 

Ae 

AE 

Distance  to  Ea 

Be 

BE 

Distance  to  E^ 

A 

A 

Distance  to  Q. 

cl 

B 

B 

Distance  to 

NOL 

NOL 

Number  of  distributed  loads  (25) 

3 

Numbers  following  explanation  of  terms  NOL--NOK  indicate 
the  arbitarily  chosen  maximum  value  of  each  for  this  program. 
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NOF 

NOF 

Number  of  point  forces  (100) 

NOR 

NOR 

Number  of  point  reactions  (10) 

NOM 

NOM 

Number  of  point  morasnts  (50) 

NORM 

NORM 

Number  of  point  reactive 
moments  (10) 

NOI 

NOI 

Number  of  bending  compliance 
elements  (25) 

NOK 

NOK 

Number  of  shear  compliance 
elements  (25) 

PDEFT 

PDEFT 

Predetermined  deflection 
at  a  reaotion 

PSLP 

PSLP 

Predetermined  slope  at  a 
reactive  moment 

NOP 

NOP 

Number  of  points  at  which 
the  flexural  quantities  will 
be  calculated. 

32 


APPENDIX  C 


General  Information  Concerning  Use  of  Program 

The  composite  program,  BEAM3  plus  subroutines,  was 
written  in  FORTRAN-60,  compiled  and  tested  on  a  Control 
Data  Corporation  l6o4  digital  computer.  It  consists  of  412 
FORTRAN  statements  requiring  about  580  cards.  The  computer 
storage  requirement  is  about  12,500  cells,  however,  this 
number  can  be  reduced  by  changing  the  dimension  statements 
which  would  of  course,  reduce  either  the  maximum  number  of 
external  loads,  elements  of  shear  compliance  or  bending  compli¬ 
ance,  or  number  of  points  at  which  flexural  quantities  could 
be  calculated. 

The  FORTRAN  statements  used  in  this  program  are:  GO  TO, 
computed  GO  TO,  IF,  STOP,  DO,  CONTINUE,  FORMAT,  READ,  PRINT, 

WRITE  OUTPUT  TAPE,  FUNCTION,  SUBROUTINE,  RETURN,  CALL,  DIMENSION, 
COMMON,  and  END. 

There  are  no  conversion  constants  incorporated  in  this 
program.  As  such  all  input  must  be  in  a  oonsistant  set  of 
units. 

"FORMAT"  for  Input  Data 

The  input  data  for  a  problem  is  read  by  the  program 
with  a  FORMAT  of  14  for  all  "fixed  point"  variables  and  F20.0 
for  all  floating  joint  variables.  The  arrangement  of  the 
data  on  the  input  oards  is  shown  as  follows: 


33 


Group 

no. 

No.  of  cards 
in  group 

Fields 

Entries 

1 

1 

F20.0 

G 

2 

1 

814 

NOP,  NOL,  NOF, 

3 

NOL 

4F20.C 

NOM,  NORM,  NOI 
A,  B,  A&>  A^ 

4 

NOF 

2F20 . 0 

Xf,  Fy 

5 

NOR 

3F20.0 

xr,  0.0,  PDEFT 

6 

NOM 

2F20.0 

V  M2 

7 

NORM 

3F20.0 

xm,  0.0,  PSLP 

8 

NOI 

4F20 . 0 

Ae*  Be»  Ea»  ^ 

9 

NOK 

4F20.0 

Ag»  Bg*  &a*  °b 

As  an  example  the  Input  data  cards  for  test  problem  #1 
In  Appendix  D  was  as  follows; 


CARD 

Is 

100 

CARD 

2; 

21  1  0 

2  0 

2  10 

CARD 

32 

0. 

100 

o 

o 

o 

I 

o 

o 

o 

i 

CARD 

4s 

0. 

0. 

.0 

CARD 

5? 

100 

0. 

oO 

CARD 

6s 

0. 

0. 

0. 

CARD 

7i 

100 

0. 

0. 

CARD 

8s 

0. 

100. 

.000000001  .000000001 
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APPENDIX  D 


Test  problems  and  solutions 

Test  problem  1 

This  problem  consists  of  a  beam  of  uniform  cross-section, 
"built-in"  at  both  ends,  and  loaded  with  a  uniformly  dis¬ 
tributed  weight  of  500#/in.  (See  Fig.  10) 


X 

X 


N. 

\1 

X 

X 


1 


Q  s  500#/ in. 


.1—1 . 1- 


100  in". 


X 

X 

X 

X 

X 

X. 


1 

El* 

K 

AG 


10~9  #in.2 


«  0 


Fig.  10 


Test  problem  2 

This  problem  is  identical  to  problem  1  except  that  the 
slope  at  either  end  of  the  beam  is  forced  to  values  other 
than  zero. 


0 

left  end 


-.002 


right  end 


.002 


35 


SOLUTION  TO  TEST  PROBLEM 


2  0<MOa>rOCTOC\JOCrroa>OCMO^f>00NOC\lO 

o  otno'-ro— olt>o 

*— •  oouir^-NicsiLmnor^cxjh-OLnLncvjpoh-uioo 

H-  OOIs-ON1J’N(\lO»-O,“O(NJNJ-M'0NOO 

O  Or-OOODrOC\ICOCOO>OCM>OOCOCOCNj^)OOOOr-0 

uj  romroo  j-o 

-J  o  O'~roun^-C7'ocsjc\jrocvjc\;oo'rwunf>0'“00 

U.  OOOOOOO—  —  —  —  r— »-0000000 

o  i  m  1 1 1 1 1 1 1 1  r  1 1  r  r  r  r  r 9 


oino^omoinoinoi/io^oinoino^o 
ocN!or-ocNJor-ocNjoc\jor-oc\jor-oc\!0 
UJ  G»“0c00'00l>00^0'-0^0>00000'-0 
a.  oooo»-oooj3,omopooj,ooof"ocoo 
o  ohois-oa'incQoouooa)ino'ONO(^o 
-J  o*—roro^fm^c\!c\j*— o^-osic^mroj'roKi^-o 
<✓>  OOOOOOOOOOOOOOOOOOOOO 
OOOOOOOOOOOOOOOOOOOOO 

I  i  I  I  I  I  I  I  I  I 


►—  ^r^r^i^f^f^ror<)r<)rofy^rororof>orON-is-r^fs-*N- 

z  o— o—ocoroco  ro  corn  co**")  corn  coo— <>*—o 

uj  o^o0'0omoi>ooroomoroo>oa''0c"0 

a:  co  cvj-"Oc\i  oo  cNjr>oc\j  oo  eg  Is- f~r--so 

o  f-O'a'O'f-inoincooooooinom^a'O'C'f- 

5:  Jd-cg—  |  I  f-»~ C\l(NC\!f-f— r-  |  I  r-CMJT 

II  I  III 


OOOOOOOOOOOOOOOOOOOOO 

o:  •  •  i  •  •  i  t  •  •  •  •  . . 

<  oooooooooo I oooooooooo 

UJ  oooooooooo  oooooooooo 

X  omooooounom  tnoino^onomo 

uo  uncgor-uncoo^LncNj  cMLnr^oc\iu^h-oc\jun 

CMCMCVi  —  r—  f—  r-  |  |  |  f— r-«— r— CMfMCM 

I  I  I  I  I  I  I 


o 

< 

o 


o 

z 

<t 


</> 


OOOOOOOOOOOOOOOOOOOOO 

OOOOOOOOOOOOOOOOOOOOO 

OOOOOOOOOOOOOOOOOOOOO 

OOOOOOOOOOOOOOOOOOOOO 

m  in  in  un  m  m  in  in  m  m  m  to  m  in  m  in  m  in  m  m  in 
I  l  I  i  l  I  I  l  l  l  II  I  i  I  I  i  l  l  I  I 


OOOOOOOOOOOOOOOOOOOOO 

OOOOOOOOOOOOOOOOOOOOO 

#•#•••••••••••••••••• 

unoinoinoLnoinounoinomOinoino 
»-  —  cm  cm  m  vj  &  j*  in  o  -o  r-  r-  co  co  o*  o*  o 


•  • 

h-N- 
•  •  *o  *o 
00-00 
OO'O'O 
OO  — 


mm.a'JTOo 
CsJCM  I  OO 

oo 
OO 
II  II  OO 

— oo 

II  II— CNOO 
—  —  OO 


f— CSjw^,  •  « 

hjrsl 

II  || 

>><<— <\l 

occtccacuu 


o 
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SOLUTION  TO  TEST  PROBLEM 


OC\IOOr0CT*O(\JOC7'f'0C'O<\JO(7'C0C7'O(\IO 
O  oi/)0^ror*omor-ro*-o</io^-roF-oino 
*“•  ooinr^-rorgminoh-ooN-omLncNj^r^unoo 
oor^oroj'h-csjo— o*— ocMh*j»m>or^oo 
CJ  OCsJCOmrorwCOrOO*— CNir-o^OCON-rnroCDCNJO 

uj  o  J’^O'LOo^oroooN^of^oornroocna'Jt^o 

-J  o*— rouico*— roLnoh-oof^oin^o^-ootr)^»-o 

u*  OOOOO^f-*-  r-^-r-  r-»—  >-^00000 

o  I  I  I  I  II  I  I  I  I  II  I  I  I  I  I  II 


o^c3inoinounou)otf)ouioi/)oinuino 
,  ocvjof^ocsjoh-ocvjocNjor^orgoh-ocNjo 

UJ  Of“UCQO'0U',f)O^-O',“O^UOOCUO'“O 
Q.  OCDOi-OOO^OfOOfOO^OOO^-OCDO 
u  otn'0^-CNic^^-3>-^CMo^  ^-=rrocrrg^oir>o 

-J  cvj^js-Lnu^  Jtfocvj^o^-Cvjro^f  js-intn  j-rocsj 
CO  OOO0OOOOOOOOOOOOOOO0O 
OOOOOOOOOOOOOOOOOOOOO 

II  II  I  I  I  I  I  I 


*-  h-h-r^h-rorofo^orororororo^rororoh-h-r^r1- 

<:  '0»-Or-r<icoro<X)rocoroaDmcX)f,ncOfOr-'0^-'0 

uj  *oo'<)0'»ynorooN,)Ofoo«oofoofoa'oa"0 

j:  'Or^^-N-rO(NCOCSjrocsjtt)c\|Ni(\jcocsjrof^r-rw<> 

O  h-imnmcMO^^  j*c\jO'^0'csiinu">inh- 

ac  roc\i*~  i  ^-r-rvjcM(\jc\j(\j^-^-  i  •— c\j**> 

III  III 


O' 

< 

UJ 

X 

00 


OOOOOOOOOOOOOOOOOOOOO 

•  •••••••••••••••••••• 


OOOOOOOOOO | oooooooooo 
oooooooooo  oooooooooo 

OLOOLOo^OLOOun  inoinomomoino 
infsjoh*mcsjo^in<\j  csKnr-ocvju^r-ocgm 
CsJCSJCM*— r— 1-r-  I  |  I  r-r-r-*~CSjCM<\J 

I  I  I  I  I  I  I 


OOOOOOOOOOOOOOOOOOOOO 
O  OOOOOOOOOOOOOOOOOOOOO 

<  «••••••••••••#••••••• 

O  OOOOOOOOOOOOOOOOOOOOO 

•J  OOOOOOOOOOOOOOOOOOOOO 

m  in  m  in  m  m  in  m  m  m  m  in  in  in  in  m  m  m  m  in  i/> 

I  I  l  I  i  I  l  I  l  i  I  I  I  I  II  i  I  I  I  l 


UJ 

O  OOOOOOOOOOOOOOOOOOOOO 

Z  OOOOOOOOOOOOOOOOOOOOO 

«1  •  •••••«•••••••••••••• 

k  moinoinoinoinoinoinoinomomo 

00  '  '"•<\IC'Jf,0*>«3^inin'0'0N’N*0DC00>0sO 
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^N- 
•  #00 
oooo 
oooo 

OONN 
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OCJ 
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II  II  *-CNOO 
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f— <\jww  •  • 

fs/fsj  I 

— I  Z  II  II 
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Test  problem  3 

This  problem  consists  of  a  beam  of  uniform  oross-section, 
"built-in1'  at  one  end  and  simply  supported  at  the  other,  and 
loaded  with  a  uniformly  distributed  weight.  (See  Fig.  11) 


V 


znoiErra 


Q,  »  l0#/ino 


40  in, 


1 

EX 

K 

AG 


,1332x10 


-7  1 


#in.2 


Fig.  11 


Test  problem  4 

This  problem  is  identical  to  problem  number  3  except 
shear  deflection  is  included  in  the  solution. 

~~  -  .159  x  10“6  l/#in. 
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SOLUTION  TO  TEST  PROBLEM 


z  Of'OCNiPOcOOCM*— Ojd‘0'0'0  0'00  C\iO 

O  OCT'OCDOin*—  Or'OinOrOrOCOOOunOCXDLO'OO 

—  OCOC^O— CMJ*  —  o*—  Or—  tnOCOrOOOCMO 

H  O^I^OO'rOr-OOOr-rO.=f  C0r-CT'r-OC7'ino 

o  oo*— roinoooroinof^coooN-o^fcvjO'Oroo 

OJ  OOOOOO^^r-rr-r-r-f-r-r-f-r-OOO 
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SOLUTION  TO  TEST  PROBLEM  4 
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Test  problem  5 

This  problem  consists  of  a  homogeneous  circular  shaft, 
three  inches  in  diameter  at  the  center  and  tapered  to  one 
and  a  half  inches  in  diameter  at  both  ends.  The  total  weight 
of  the  shaft  is  100#  and  a  uniformly  distributed  load  total¬ 
ing  200#  1 8  placed  on  the  central  section.  The  shaft  is 
supported  by  a  uniformly  distributed  reaction  at  each  end. 
Points  A  and  B  are  points  of  zero  deflection.  This  problem 
was  approximated  for  the  program  as  shown  in  Fig.  12. 


Fig.  12 


Weight  of  shaft: 


0"  -  8"  -0.5315#/in. 

8"  -  24"  -0.5315  to  -2.081#/in.  (linear) 

24"  -  36"  -2.081#/in. 


Bending 

compliance  (—)  ; 

0"  -  8" 

0.1420x10  ^ 

l/#-in?  , 

8"  -  12" 

0. 1420x10" ? 

to 

0.0583x10"? 

l/#-l nf 

( linear) 

12"  -  16" 

0.0583x10-9 

to 

0.0281xl0”y 

l/#-ln| 

(linear) 

16"  -  20" 

0.0281x10"? 

to 

0.0152x10"? 

l/#-ln| 

(linear) 

20"  -  24" 

0.0152x10"? 

to 

0. 0089x10" b 

(linear) 

24"  -  36" 

0.0089x10'  5 

l/#-in? 

Shear  compliance  {]&)  -  0. 
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SOLUTION  TO  TEST  PROBLEM 
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Teat  problem  6 

This  hyperstatic  problem  consists  of  a  beam  with  a  very 
unlikely  configuration  as  shown  in  Fig.  13 .  The  bending  com¬ 
pliance  varies  as  shown  in  Fig.  14.  Shear  deflection  has  been 

K 

oml tted,  1 . e . 9  —  0  s 


Fig.  14 


Test  problem  7 

This  problem  is  identical  to  problem  number  6  except  the 
beam  was  turned  end  for  end. 
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SOLUTION- TO  TEST  PROBLEM  6 
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SOLUTION  TO  TEST  PROBLEM  7 
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APPENDIX  E 


Subroutine  GAUSS2 

This  routine  is  a  modified  form  of  the  GAUSS2  used  by 
C.  B.  Bailey  /7/  in  his  program,  LINEQN,  for  solving  simul¬ 
taneous  linear  equations  by  Gaussian  elimination.  Bailey 
was  solving  the  matrix  equation,  Ax^B,  with  a  possibility 
of  50  vectors  B  for  each  matrix  A,  but  this  program  requires 
only  one  vector  B  with  each  matrix  A.  Also  since  the  maximum 
number  of  equations  in  BEAM3  is  22,  the  size  of  matrix  A  was 
reduced  from  100x101  to  22x23. 

At  each  step  of  the  elimination,  the  value  of  the  diagonal 
element  is  compared  to  a  defined  11  zero ".  If  it  is  smaller 
than  this  quantity,  the  matrix  is  considered  singular  and 
an  error  output  is  returned  to  the  calling  program.  This 
is  the  "matrix  singular"  test  used  by  BEAM 3 . 


APPENDIX  F 


Listing  of  program  and  subroutines 


Name 

Page 

BEAM3 

48 

LOAD 

54 

SHEAR 

55 

MOMENT 

56 

SLOPE 

58 

DEFLECT 

61 

FTON  and  FT1N 

64 

ENTON  and  ENT1N 

65 

DITON  and  DITIN 

66 

GAUSS 2 

67 
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PROGRAM  BEAM3 

DIMENSION  W(  1001 )  , V ( 1 001  ) , AM ( 100 1 ) , Y 1 ( 1 00 1 ) , Y { 1 00 1  ) , XX(  100 1  )  , 

1  A(25) , B( 25) ,QA(25) ,QB(25) ,XFO( 100), FY( 100) ,XR( 1 0),RY{ 10) , 
2XM(50)  ,AMZ  (50)  ,RXMMO)  ,RAMZ(  1 0 ) , AE ( 25 ) , BE ( 25 ) , E A { 25 ) , EB ( 2  5 ) , 
3AA(22,23) , XA(22) ,PDEFT(1  0),PSLP(10) 

4  , 1  TITLE ( 10)  ,GA (25  )  ,GB(25) , AG (25) ,BG (25 ) 

COMMON  A,B ,QA,QB, XFO, FY, XR , RY , XM , AMZ , RXM , RAM Z , AE , BE, EA,E3,N0L, 
1  NOR, NOF,NOM, NORM, NOI , NO K , G A , GB , AG , BG 
READ  70, C 

70  FORMAT  (F20.0) 

READ  71 , NOP , NOL , NOF , NOR, NOM , NORM, NOI.NOK 

71  FORMAT (814) 

IF(NOL)  51,51,50 

50  READ  72 , ( A ( I  ),B(I),QA(I)  , QB ( I ) ,1  =  1, NOL) 

72  FORMAT ( 4F20 . 0 ) 

51  IF( NOF) 53, 53,52 

52  READ  73,  (  XFO(N)  ,FY(N)  ,'4=1, NOF) 

73  FORMAT  (2F20.0) 

53  IF( NOR) 55, 55,54 

54  READ  74, (XR(M),RY(M) ,PDEFT(M),M=1 ,NOR) 

74  FORMAT  (3F20.0) 

55  I F ( NOM ) 57 ,57 ,56 

56  READ  73,  (XM( K) ,AMZ(K) ,K=1 ,NOM) 

57  I F ( NORM ) 59,59,58 

58  READ  74, (RXM(KA) , RAMZ ( KA ) , PSLP( KA) ,KA=1,N0RM) 

59  IF(NOI ) 6 1 ,61,60 

60  READ  72, (AE(  IA), BE ( I  A) ,E A( IA) ,EB( IA) , I A= 1 , NO  I ) 

61  I  F  (  NOK ) 1 , 1 , 62 

62  READ  72,  < AG( IB) ,BG( IB),  GA(IB) ,GB( IB), IB=1,N0K) 

1  NORl  =  NOR+1 

C  1  =  0. 

C2=  0. 

NORM 2  =NOR  M  +  2 
N0RM3  =NORM  +  3 
NOQ  =  NOR  +NORM 
NCQ1  =  NOQ  +  1 
N0Q2  =  NOQ  +  2 
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N0Q3  =  NOQ  +  3 
EPSX= 1 0 . ** ( -8 ) 

DO  2  K=  1  ,  NCR 

2  AA(1,K)  *  1. 

DO  3  K=  NCR  1, NOQ 2 

3  AA(  1  ,K  )  =  C. 

AAA=G. 

IF (NOF) AS ,45,35 
35  DO  4  N  = 1 , NC F 

4  AAA  =  AAA  +  FY ( N) 

45  I F ( NOL ) 48,48,46 

46  DO  5  1=1, NCL 

5  AAA  =  AAA  +  ( ( QA (  I ) +QB ( I  1) /2 . ) *  (B(I)-A(I)) 

48  AA ( l,NOQ3)=-AAA 

DO  6  K  =1  ,  NCR 

6  AA ( 2  , K )  =  C-XR(K) 

IF{NCRM)9,9,7 

7  DO  8  K=NOR 1 , NOQ 

8  AA { 2 , K  )  =  1  . 

9  DO  1C  K=NCQ 1 , NOQ2 

10  AA ( 2  ,K  )  =  0. 

CALL  MOMENT  ( C ,A AA ) 

AA ( 2 ,N0Q3 )  =-AAA 
IF  (  NORM  117,17,11 

11  DO  16  J=3, N0RM2 
DO  13  K= 1 , NCR 
AA ( J , K ) =0. 

DO  12  I A  = 1 »NOI 

AAA  =  EA(  I  A ) *ENT  0N(  RXM ( J-2 ),AE(IA),XR(K),1)  -  EB(IA) *ENTON ( RXM ( J-2 )  , 
1 BE( IA) ,XR(K), 1 )  +((EB(IA)-EA(IA))/(Bc(  I  A )-AE ( I  A)  )  )  * ( ENT  1 N ( RXM ( J-2 ) 
2, AE ( I  A) ,XR(K) , 1 )- ENT 1N( RXM( J-2) , B  E {  I A  )  ,XR(K),  1 ) ) 

12  A A ( J , K )  =  A A ( J , K 1  +  AAA 

IF  (NCK 1  13,13,900 

900  DO  91 C  I B= 1 , NOK 

IF  (  NORM-1  1902 ,901  ,9C2 
902  IF( J-N0RM21901  ,904,901 
9 C 1  X  =  RXM ( J-2 )  -EPSX 
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AAA  =  GA (  IB  )  *  FTO  N (  X  , AG (  I B  )  , XR ( K ) , 0 )  -  GB(IB)*  FTOMX 

1  , BG( 18), XR(K)  ,0)  +  {  (GB{ IB)-GA(IB) )/{BG(IB)-AG( IB)  ) )*  (  FT  IN  ( 

2  X  , AG(IB) ,XR(K) ,0)-  FT)  N  (  X  ,  BG  (  IB  )  ,  XR(  K  )  ,  0  )  ) 

GO  TO  910 

* 

904  X=RXM ( J-2 )  +EPSX 

AAA  =  G  A (  16)*  FTO  N (  X  , AG ( I B )  , XR { K ) , 0 )  -  GB ( I B  )  *  FTCN(X 

1  , EG  C IB) ,XR(K)  ,0)  +  (  (GB( IB)-GA(IB))/(BG(IB)-AG( IB) ))«  (  FT  1 N ( 

2  X  , AG(IB) ,XR(K) ,0)-  FT  1 N {  X  , BG ( I B ) , XR { K ) , 0 ) ) 

910  A  A  (  J  ,  K' )  =  A  A  (  J  ,  K  )  -AAA 

13  CONTINUE 

IF( N0C-N0R1  )  1  55,  1  35,  135 
135  DO  15  K=N0R1,NCQ 
AA{J,K)=  0. 

DO  14  I  A  =  1 ,NOI 

AAA=  EA ( IA) *ENTON (RXM( J-2)  ,AE(  I  A ) ,RXM (K-NOR )  ,  C )  -EB(IA)*  ENTCN 
1  ( RXM (J-2) , BE(IA)  ,RXM(K-NOR) , 0)  +  ( ( EB ( I  A ) -E A ( I  A ) ) / ( BE ( I  A )-A E ( I A ) ) 
2)  *(ENT1N(RXM( J-2 ) ,AE( IA) ,  RXM ( K-NOR ) ,0)  -  ENT  1 N { RX M { J-2 ) , B E (  I A ) , 

3  RXM { K-NCR ) ,0) ) 

14  AA{ J,K)  =  AA( J,K )  +  AAA 

15  CONTINUE 

155  AA ( J , NOG  1  )«1. 

AA ( J , N0C2 ) =  0. 

-  I F ( J-N0RM2)  156,  157, 156 

156  X=RXM ( J-2 )  -  EPSX 
CALL  SLOP  E ( X , AAA) 

GO  TO  16 

157  X=RXM ( J-2 )  +  EPSX 
CALL  SLOPE! X,AAA) 

16  A A( J , NOG  3 )  =- AAA  +  PSLPCJ-2) 

17  00  2C  J=NORM3,NOC2 
DO  IE  K= 1 , NOR 

AA ( J , K )  =0. 

DO  181  I  A=  1  , NO  I 

AAA  =  EA(IA)*  DI TON( XR( J -NORM 2 ),AE(IA),XR(K),  1)  -EB ( I  A) *CITCN 

1  ( XR ( J-N0RM2 ) » BE(  I A  )  , XR ( K ) , 1 )  + ( ( E B ( I  A )-EA ( I A ) ) /  ( BE (  I A ) -A E (  I A ) )  ) 

2  *(DIT1N(XR( J-N0RM2), AE( IA) ,XR(K) , 1 )  -  D I T 1 N ( XR { J-N0RM2 ) , 8 E ( I A ) , 

3  XR( K  )  , 1 ) ) 
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181  A A ( J  *  K )  =  AA(J,K)  +  AAA 
IF  { NOK ) 1 8, 18,920 
920  DO  930  18=1 ,NOK 

AAA  =GA (IB)*  ENTON{ XR ( J-N0RM2 ) » AG(  IB) ,  XR ( K) ,  0)  -  G8 ( IB ) *E  NTON 

1  (XR{ J-N0RM2 ), BG( 18)  ,XR( K) ,0)  +( (GB( ID)-GA( IB) )/(BG(  I B ) —  A G ( IB)  ) ) 

2  * ( ENT1 N ( XR ( J-N0RM2 ) , AG(  IB)  ,XR(K),0>  -ENT1N t XR( J-N0RM2) , BG( IB) , 

3  XR( K) ,  0)  ) 

930  AA( J ,K )  =AA { J, K )  -AAA 

18  CONTINUE 

I F (NOO-NOR 1  )  195,185,  185 
185  DO  19  K=NOR 1 ,NOQ 
AA ( J , K )  =  0. 

DO  19  I A= 1 , NOI 

AAA  =  EA ( I  A )*  DITONIXRC J-N0RM2)  , AE ( I  A) , RXM( K-NOR ) , 0 )  -  E  B { I  A ) * 

1  DITONt  XR { J-N0RM2 ) , BE( I  A ) , RXM { K-NCR )  ,0 )  + ( ( E B C I  A ) -E A ( I  A ) )  /  { BE ( I A 

2  )-  AE (  I A ) ) )  *  ( DI T 1 N ( XR ( J-N0RM2 ) , AE ( I  A) ,RXM( K-NOR) ,0)  -  DITIN 

3  (XR ( J- NOR M2 ) , BE ( I  A)  , RXM ( K-NOR ), 0 ) ) 

19  AA ( J, K )  =  A A ( J , K )  +  AAA 
195  A A { J , N0Q1 )  =XR( J-N0RM2) 

AAC  J , N0Q2 )  =  1. 

CALL  DEFLECT  { XR  C J-N0RM2) , AAA) 

20  AA(J,N0Q3)  =-AAA+  PDEFT ( J-N0RM2 ) 

-  CALL  GAUSS2 (N0Q2, .00000001  ,AA,  XA,K1 ) 

GO  TO  (22,2 1 ) ,K  1 

21  WRITE  OUTPUT  TAPE  4,77 

77  FORMAT  ( 16H  MATRIX  SINGULAR) 

STOP 

22  DO  23  K= 1 , NOR 

23  RY  (  K)  =  X A(  K ) 

I F ( NORM ) 26,26,24 

24  DO  25  K=NO Rl ,NOQ 

25  RAMZ ( K-NOR )  =  XA(K) 

26  Cl  =  XAtNOQl ) 

C2  =  XA ( N0Q2) 

X=0. 

NOPE  =N0P-1 
FNOP  =  NOP- 1 
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DO  27  J=1,NCPE 
CALL  LOAD  (X,W(J ) ) 

CALL  SHEAR  (X,V( J ) ) 

CALL  MOMENT  ( X, AM( J )  ) 

XI  =X+EPSX 

CALL  SLOPE (  X  1  ,  Y  1  ( J) ) 

Y  1  {  J ) = Y 1 ( J  )  +C1 
CALL  DEFLECT  (X,  Y(  J)  ) 

Y(J)=Y(J)+C1*X+C2 
XX ( J )  =  X 
27  X  =  C/FNOP  +  X 
XX (NOP  )=X 
'  X=X-3.*EPSX 

CALL  LOAD (  X  ,  W  ( NOP  )  ) 

CALL  SHEAR(X, V(NOP) ) 

CALL  MOMENT  (X»AM( NOP)  ) 

CALL  SLC° E (  X ,  Y 1 (NOP) ) 

Y 1 ( NO P  )  =Y  1 ( NOP  )  +C  1 
CALL  CEFL  ECT ( X, Y( NOP)  ) 

Y(NOP)=Y(NOP)+Cl*X  +C2 
WRITE  CUTPUT  TAPE  4,76 

76  FORMAT  (  3X , 1 H J ,7 X , 4HL0AD, 1  OX, 5HSHEAR , 1  OX , 6HMCMENT , 8X , 5HSLCP E , 

-  1  9X, 1CHDEFLECTION ,8X, 1 HX ) 

WRITE  OUTPUT  TAPE  4 , 7 8 , ( J , W ( J ) , V ( J ) , AM ( J )  , Y 1 ( J ) , Y( J ) , XX ( J )  , J  =  1  , 
1  NOP  ) 

78  FORMAT! 14, 6E1 5.8) 

IF  (NCR  ) 2 82 , 282,280 

280  PRINT  281  , (M,RY(M ) , M= 1 ,NOR) 

281  FORMAT  (  3HRY ( , 12 , 2H ) = , E20. 8 ) 

282  IF(N0RM)288,288,283 

283  PRINT  284, (KA,RAMZ(KA  ),KA  =  1 , NORM) 

284  FORMAT  (5HRAMZ ( , I  2, 2H )  =  , £20. 8) 

288  PR  INT  286, C  1 

286  FORMAT (5X,3HC1*,E20.8) 

PRINT  267, C2 

287  FORMAT (5X,3HC2=,E20.8) 

285  CONTINUE 


STOP 

END 


* 


-*-■ 


i 
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SUBROUTINE  lOAO(X,W) 

DIMENSION  A (25) ,B (25) ,QA{ 25) ,CB(25) ,XFC( 100) ,FY( 100) , XR( 10) ,  RY( 1C) 
1 »XM( 50  f AMZ(50) ,RXM( 10) ,RAMZ( 10),AE(25),RE(25),EA(25),EB(25) 
2,GA(25) ,GB(25),AG( 25),BG(25) 

COMMON  A, B,QA,QB, XFO.FY, XR  ,RY  ,XM, AMZ, RXM, R AMZ , AE , B E »  E A , E B t NOL , 

1  NOR»  NOF, NOM ,NORM  »  NCI »  NOK  » GA ♦ GB  » AG , BG 
EPS  X= 10 .** (-8 ) 

w=o .  . 

1  00  5  I-  1 1  NOL 

2  IF(X+EPSX-B(I))3, 3,5 

3  IF(  X+EPSX-A(  I  )  )5»  4,Ij 

4  WW=QA ( I ) + ( ( QB( I )-QA ( I ) )*(X-A(I  ) ))/(B(I)-A(I)) 

w=w+ww 

5  CONTINUE 
RETURN 
END 
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SUBROUTINE  SHE  AR( X ,  V ) 

DIMENSION  A (25 ) ,B (25) ,QA(25) ,QB(25) ,XFC(  ICC) , FY (100) , XR ( 1C) , RY (1C ) 
I*XM(5C)  ,  AMZ ( 50 ) ,R  XM  ( 1 C ) »  RAMZ ( 10  )  ,  AE(  2 5  )  ,  BE  (25  ) ,  E  A(  25  )  ,  EB  (  2  5 ) 
2,GA(25),GB(25),AG(25)»BG(25) 

COMMON  A, B,QA,QB, XFO, FY , XR, RY , XM , AMZ , RXM, R AMZ , AE ,B E , E A , E B , NOl , 

1  NOR , NOF , NOM , NORM ,NCI,NOK,GA,GB,AG»BG 
EPSX=10.**(-8) 

V=0  . 

IF ( NOL ) 7, 7 , 1 

1  DO  6  1=1, NOL 

IF ( X-A( I ) ) 6,2,2 

2  I F  (  X+EPSX-B(  I ) )4, 4,3 

3  VV=  QA( I ) * ( X- A ( I )  )  +  ( (GB( I )-QA( I ) ) * ( ( X- A ( I ) ) ** 2 ) ) /( 2 . « ( B ( I)-A  ( I  )  )  ) 
l-QBl I ) » ( X-B ( I ) )  - ( (Q8( I )-QA ( I ) ) *( ( X-B( I ) ) **2) )/ (2.  *(B( I >-A( I ) ) ) 

GO  TO  5 

l4  VV=+Q  A  (  I  )  *(X-A(  I)  )+(  (GB(  I  )-QA(I  JO*  (  (X-A(  I))»*2))/(2.*(B(I)-A(I))) 

5  V=V+VV 

6  CONTINUE 

7  IF(N0F)10, 10, 75 
75  DO  9  N=  1  ,  NOF 

IF ( X+EPSX— XFO ( N) )  9,9,8 

8  VV=FY(N) 

v=v+vv  ' 

9  CONTINUE 

10  IF (NOR  113,13,105 
105  DO  12  M= 1 ,NOR 

I F ( X+  EPSX-XR ( M ) 112,12,11 

11  VV=  R  V  (  M  ) 

V=V+VV 

12  CONTINUE 

13  CONTINUE 

RETURN 

END 


i 


55 


SUBROUTINE  FOMENT ( X, AM ) 

DIMENSION  A(25),B( 25) ,QA(25) ,QB(25) ,XFG( 100) ,FY(100)  ,XR(  1  C),  RY(IO) 
1 «  XM ( 5C ) ,  AMZ ( 50 } »R  XM ( 1 0 ) »  R  AMZ ( 10),AE(25),BE(25),EA(25),EB(25) 
2,GA(25),GB( 25) ,  AG (25) »BG(25) 

COMMON  A, B,QA,QB, XFO , FY , XR , RY , XM, AMZ, RXM , R AMZ , AE , B E , EA , EB , MOL , 

1  NOR, NGF , NOM ,NORM , NOI »  NOK , GA, GB»  AG  »  8G 
EPSX=10.«*(-8 ) 

AM  =  0. 

IF(NOL)8,8,2 

2  00  7  1=1 , NOL 

IFIX-Al I ) )7, 3 , 3 

3  IF  (  X+EPSX-B(I  )  )5,  5,4 

4  AMM  =  ( +CA(I)/2. )*(  ( X-A( I) ) **2 )  +  ( (QB( I ) -QA  C I )  ) * ( (X-A ( I )> **3 ))/ ( 6 .« ( 
1BII )-A(I) ) )-(OB(I )/2. )*( ( X-B( I ) 1**2 >-  (QB(  I  >-CA ( I )  )  * ( ( X-B (I)  )**3> 
2/  (  6. * ( B ( I ) —A (  I)  )) 

GO  TO  6 

5  AMM=( +  Q  A  (  I )/2. )*(  ( X — A (  I) )**2)+( ( QB ( I ) -QA ( I } }* ( (X-A ( I ) )**3) )/ (6. * 

1 (B(  I )  —  A ( I  )  )  ) 

6  AM= AM+AMM 

7  CONTINUE 

8  IF(NOF) 12, 12,9 

9  DO  11  N=1 , NOF 

-  IFC X+EPSX-XFOIN))  11,11,10 

1 

10  AMM=  F Y ( N ) * ( X-XFO (N) ) 

AM  =  AM  +  AMM 

11  CONTINUE 

12  IF(NOR) 16, 16,13 

13  DO  15  M= 1 , NOR 
IF(X+EPSX-XR(M) >15,15,14 

14  AMM=  R Y ( M ) * ( X— XR( M)  ) 

AM=AM+AMM 

15  CONTINUE 

16  IF (NOM 128 ,28, 17 

17  DO  19  K= 1 , NOM 

I  F  (  X  +  EPSX-XM(  K )  119,19,18 

18  AMM= AMZ ( K 1 
AM= AM+AMM 
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19  CONTINUE 

28  IFtNORV)  32,32,29 

29  CO  31  KA= 1 ,NORM 

I F  C X+EPSX-RXMIKA)  )31,31,30 

30  AMM  =  R AMZ ( KA ) 

AM= AM+AMN 

31  CONTINUE 

32  CONTINUE 
RETURN 
ENO 


\ 
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SUBROUTINE  SLOPE!  X,Y1  ) 

DIMENSION  A(25 ) ,8 (25) ,QA(25) , QB  C  2  5 ) ,XFC(  100) ,FY(100) , XR( 10) ,RY(  10  ) 

1 , XM(5C) ,AMZ(50) ,RXM( 1 0 ) , RAMZ ( 1 0 ) , AE (25 ) , B E ( 25 ) , E A( 25 ) , EB { 2 5 ) 

2, GA( 25) ,GB(25), AS (25) ,BG(25) 

COMMON  A, B,QA, QB, XFO, FY, XR , RY, XM, AMZ , RXM , RAMZ , AE ,B E, EA*  EB» NOL  , 

1  NOR,NOF,NON, NORM  » NO  I , NOK , GA , GB , AG, BG 
Y 1  =  0  . 

00  17  I  A=  1 , NO  I 
IF(X — AE(  IA)  )  17, 17 ,2 

2  IF(N0L)5,5,3 

3  00  4  1=1, NOL 

20  EYY  =  EA( IA)M ( QA ( I ) /2. ) *ENTON( X, AE( I  A) ,A( I), 2)  + (  { QB ( I ) -Q A (  I  )  )  / 

1  (6.  *  (  e  (  I  )-A(  I  )  )  ) )  MENTON  (X,AE(  I  A)  ,  A(  I  ) ,  3  J-ENTON  (  X,  AE  (  I A  )  ,B  ( I  )  ,  3  \  ) 

2  -{QB (I) /2. )*ENTON( X, AE( I  A) ,B(  I) ,2 ) )  +( ( EB ( I A)-EA (IA))/(BE(IA)- 

3  AE ( I A ) ) ) * ( (0A(I)/2.)»ENT1N(X,AE(IA),A(I),2)  M  ( GB  ( I ) -CA  (  I  )  )  / 

4 ( 6 . * ( E ( I ) — A ( I ) ) ))* (ENT1N( X, AE( I  A) , A ( I ) ,3 )-ENT INC  X,  AE ( I  A ) , B ( I  ) ,3)  ) 

5  -(QB(  I  )/2.  )*ENT1N(X,AE(  IA),B(I)  ,2)  ) 

Y1=Y1 +EYY 

21  EYY=  -EB (  IA)*((QA(I)/2.)*ENT0N(X,BE(IA),A(I),2)  + (  ( GB (  I  ) -C A (  I  )  )  / 
1(6.*(B(I)-A(I)))) MENTON(  X,BE(  I  A) , A ( I ) , 3 )- ENT  ON ( X, BE ( I  A ) , B ( I ) , 3 ) ) 

2  - { QB ( I )/2. ) *ENTON ( X, BE( IA),B(I),2))  -( ( EB ( I  A ) — EA ( I A ) ) / { BE { I A )- 

3  AE(IA)))*( (QA(I)/2.)*ENT1N(X,BE( IA),A(I),2)  + ( ( CB ( I ) -QA ( I ) )  / 
4'(6.*(  B(  I  )v- A  {  I  )  )  ) )  *  (  ENT  IN  (  X,  BE  (  I A )  ,  A  (  I  )  ,  3  J-ENT 1 N  (  X,  BE  (  I A  )  ,  B  ( I  )  ,  3  )  ) 

5  -(Q8(I)/2.)*ENT1N(X,BE(  IA),B(I),2)) 

Y 1=Y 1 +EYY 

4  CONTINUE 

5  IF(N0F)8,8,6 

6  DO  7  N= 1 , NOF 

22  EYY  =  E  A  {  IA)*FY(N) *ENTON(X,AE( IA),XFO(N) , 1 )  -E B ( I  A ) *FY ( N ) *ENTON ( X  , 

1 BE( I A ) , XFO ( N ) , 1 )  +  ( (EB( IA)-EA( IA) )/(BE( IA)-AE(IA) ) )*FY(N)« 
2(ENT1N(X, AE( I  A) ,XFO(N  ), 1 )  -  ENT  1 N ( X, BE ( I  A  ) , XFO (N ) , 1 ) ) 

Y1 =Y1 +EYY 

7  CONTINUE 

8  IF(NOR)  1  1  ,  11,  9 

9  00  10  M= 1 , NOR 

23  E YY=  EA(IA)*RY(M)*ENT0N(X,AE(IA),XR(M),1 ) -EB ( I  A ) *R Y ( M ) *  ENTON ( X , 

1  BE ( IA),XR(M),1)  + ( (EB(IA)-EA( IA) )/(BE(IA)-AE( IA) ))*RY(M)  * 
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2(ENT1N(X,AE(  I A  )  ,X  R  ( M  )  ,  1  )  -  ENT1N(X,BE(  I A  )  ,  XR  ( M  )  ,  1  ) ) 

Yl =Y1 +  EYY 

10  CONTINUE 

11  I F ( NOR  )  14,14,12 

12  00  13  K* 1 ,NOM 

» 

24  EYY=  E A  < I  A) *AMZ(K )*ENTON( X, AE ( I  A) ,XM( K) , 0)  -EB< I  A) *AMZ (K ) » ENTON ( X , 

1  BE  ( I A  )  ,XM(*),0)  +  (  (EB(IA)-EA(IA)  ) /(BE(  I A  )-*AE  ( I A  ) )  )•  AMZ(K)  * 

2  (ENT1N(X,AE( IA),  XM(K) ,0)-  ENT1  N t X* BE ( I A ) , XM ( K ) , 0)  ) 

Yl=Yl+EYY 

13  CONTINUE 

14  I  F ( NORM )  17,17,15 

15  DO  16  KA=  1 , NORM 

25  EYY=E A (  IA  )  *RAMZ <K A ) *ENTON ( X ,  AE  ( I A ) , RXM (KA ) , 0 ) -EB ( I  A ) * RAMZ ( K A ) * 

1  ENT0N(X,BE(IA),RXM(KA),0)  +  < < EB ( I A) -E A( I A ) ) / < BE( I A ) -AE( I A ) )  )  • 
2RAMZIKA)* (ENT  IN (X ,AE( IA) ,RXM(KA) ,0)-  ENT  1 N ( X, BE ( IA ) , RXM ( KA ) ,01) 

Yl  =Y1 +  EYY 

16  CONTINUE 

17  CONTINUE 
IF(NCK)595, 595,500 

500  DO  59C  16=1, NOK 

IF  (X  —AG (IE) 1590,590,509 

m 

509  IF  (NCL 1530, 530,5 10 

510  00  52  C  1=1, NOL 

26  EYY  =  GA( IB)*(QA(  I ) *FTON ( X, AG ( I B ) ,A(I ),1 )  ♦  (<QB(I)  -Q A ( I  )  )  / 

1  ( 2. * ( B ( I )  —  A ( I  1 ))  )*(FTON(X,AG( IB) ,A(I ) , 2 l-FTON ( X ,AG ( I B ) , B ( I ) ,2)1 

2  -Q8( I )*FTON( X»AG ( IB  ) ,B( I ) , 1 ) )•  +  ( ( GB ( IB )-GA { I B) )/ < BG( IB )-AG ( I B  )  )  ) 

3  *(QA( I )*FT1N(X,AG( IB) ,A(I) , 1 )  +  ( ( QB ( I ) -Q A (I) ) / ( 2 . * ( B ( I ) - A ( I ) ) ) ) 
4*  ( FT1N(X,AG( IB), A ( I ) ,2)  -FT  1 N( X, AG ( I B ) , B ( I ) , 2 ) )  -  QB ( I ) *  FT  1 N 

5  (X,AG(IB  ),B(I)  ,2)  ) 

Yl  =Y 1  -EYY 

27  EYY  =-GB( IB)*(QA( I ) *FTON ( X , BG ( I B ) , A( I ),1 )  +  <(QB(I)  -QA(I))  / 

1  (2.*(B(I)-A( I))) )*(  FT0N(X,BG(IB),A(I),2)-FT0N(X,BG(IB),B(I),2)) 
2-QB ( I )*  FTONI X,BG (IB),B(I),1))-  ( ( GB ( I B ) -GA ( I B ) ) /( BG( IB)-AG(  IB) ) ) 
3* (QA (  I  )  *  FT1N(X,BG(IB),A(I)»1 )+  ( ( QB( I )-QA ( I ) )/ (2. * (B ( I )-A ( I  )  )  )  ) 
4*(  FT1N(X,BG( IB), A(I ) ,2)-  FT  1 N( X, BGII B) , B ( I ) , 2) )  -  QB( I )*  FT1N 
5  (X,BG( ie  ), B( I )  ,2  )  ) 

Yl  =Y1  -EYY 
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520  CONTINUE 

530  IF(NOF)56C»560»540 

540  DO  55 C  N=1,N0F 

28  EYY=GA ( IB )  *FY ( N )*  FTON ( X , AG ( I B ) , XFO ( N ) , 0 )  -G B ( I B ) * FY ( N ) *  FTON(X, 
1BG(  18)  ,XFO(N)  ,0)  +  (  (GB,(  I  B  > -G  A  (  18)  )  /  (  BG  (  IB  ) -AG  ( I  B  )  )  )  *  FY  ( N )  * 
2 (  FT  1 N ( X, AG( IB) ,XFO(N) ,0)  -  FT1 N ( X, BG ( I B ) ,XFC <N ) , C ) ) 

Y 1  =Y1  -EYY 
550  CONTINUE 
560  IFINOP )59C,590,570 
570  DO  58 C  h‘  *  1  » N 0 R 

29  EYY  =GA (IB) *RY ( M) *  FTON ( X , AG ( I B ) , XR (M ) ,0  )  -GB ( I B ) * RY ( N ) *  FTON(X, 
1 BG ( IB)tXR(M).O)  +  ( (GB(IB)-GA(IB) )  /  ( BG ( I B ) -AG ( IB ) )  )  *  RY ( M  )  * 

2 (  FT1N(X,AG(IB),XR(M),0)  -  FT  1 N ( X, BG ( I B ) , XR ( M ) , 0)  ) 

Yl  =Y  1  -EYY 
580  CONTINUE 
590  CONTINUE 
595  CONTINUE 
RETURN 
END 
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•  SUBROUTINE  DEFLEC T ( X , Y ) 

DIMENSION  A(25) ,B (25) ,QA(25) ,QB(25),XFO( 100) ,FY( 100),XR( 10)»RY( 10) 
1,XM(50),AMZ(50) »RXK(10),RAMZ(  10  ) , AE ( 2 5) , BE (25 ) , EA( 25 ) , EB( 25 ) 
2»GA(25)»GB(25)»AG(25),BG(25) 

COMMON  A,B»QA,QB, XFO , FY , XR ,RY , XM, AMZ, RXM , R AMZ , AE ,B £ , EA , EB , NOL , 

1  NOR,NOF,NOM,NORM,NOI,NOK,GA,GB,AG,BG 
Y=0. 

DO  17  I A  =  1 » NO  I 
IF(X-AE(IA) )17,  17,2 

2  I F ( NOL 15,5,3 

3  DO  4  1=1, NCL 

20  EYY  =  EA( IA)*((QA(I)/2.)*DIT0N(X,AE(IA),A(I),2)  +(  {  QB ( I )-QA ( I ) )/ 

1 (6.*(B( I )-A( I ) ) )) *(DITON(X,AE( I  A) ,A( I ),3)-DIT0N( X,  AE  ( I  A  >  ,  B  (  I ) ,  3  )  ) 

2  -(QB(I)/2.)«DIT0N(X,AE(IA),B(I),2))  +  (  (  EB  ( I  A  1  -EA  (  I A  ))  /  (  BE  (  I A )  - 

3  AE( IA) ) )*{ (QA( I)/2.)«DIT1N(X,AE( IA),A{ I  1,2)  + ( (QB { I ) -QA ( I ) )  / 
4(6.*(B(I)-A(I))))*(DIT1N(X,AE(IA),A(I),3)-DIT1N(X,  AE(IA), 6(11,3)) 
5-(QB( I )/2. )*DIT1N(X,AE( I A),B(  I  ) ,2) ) 

Y=Y  +  E  YY 

21  EYY=  -EB(IA)*((QA(I)/2.)*DIT0N(X,BE(IA),A(I),2)  +(  ( QB ( I ) -Q A (  I )  )  / 

1 (6.*( B( I )-A( I ) ) ))*(DlTON(X,BE( I  A) , A(I  ) , 3 ) -DITON ( X,  BE ( I A ) ,  B  (  I ) , 3 ) ) 

2  -(QB ( I )/2. )*DITON(X,BE( IA),B( I ) ,21)  -((EB(IA)-EA(IA))/(BE(IA)- 

3AE(  I A  )  )  )*(  (QA(  I)/ 2.  )*0IT1N(X,BE(.IA)  ,  A  ( I  )  ,2)  +  (  (QB(  I  )-QA(  I  )  )  / 
4(6.*(E(I)-A(I))))*(DIT1N(X,8E(IA),A(I),3)“DTT1N(X, BE(IA),B(I),3)) 

5  - <  QB ( I )/2. )*DIT1N(X,BE( IA),B( I ),2) ) 

Y=Y+EYY 

4  CONTINUE 

5  IF(N0F)8,8,6 

6  DO  7  N= 1 , NOF 

22  EYY=  EA( IA)*FY(N) *DITON(X,AE( IA) ,XFO( N), 1  )  -EB ( I  A) *FY( N) *D  ITON ( X , 

1  BE( IA) ,XFO(N), 1)  +  ( ( EB ( I A ) — E  A ( IA) )/(BE( I A )  —  AE ( IA ) 1 )*FY(N)* 
2(DIT1N(X,AE(IA),XF0(N),1 )  -  DIT IN ( X ,B E ( I  A ) , XFO( N ) , 1 )) 

y=y+eyy 

7  CONTINUE 

8  IF ( NOR  111,11,9 

9  DO  10  M=1 , NOR 

23  EYY=  EA( IA)*RY(M)  *DITON( X, AE( IA) ,XR(M), 1  )-EB(  IA)*RY(M1*  DITON(X, 

1  BE ( I A ) , XR ( M ) , 1 )  +((EB(IA)  — EA( IA) )/(BE(IA) — AE (IA)))*RY(M)  * 
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2(0IT1N(X,AE( I  A  )  ,  X  R  (  M  ) , 1 )-  DITlN(X,BE(IA),XR(M),l)) 

Y=Y+EVY 

10  CONTINUE 

11  I F ( NON ) 1 4 » 1 4»  12 

12  DO  13  K= 1 , NOM 

24  EYY=EA( IA)*AMZ(K) *D I  TON ( X , AE ( I A ) , XM ( K ) , 0 )  -EB( I  A ) * AMZ ( K ) * D I TON ( X , 

1  BE ( I A ) t  XM ( K  )  t  0  )  +  ( (EB(I A)-EA(IA) )/(BE( IAJ-AE (IA))  )*AMZ(K)*  ! 

2  (DIT1NU, AE( IA),  XM{K) ,0 )-  DI T  1 N ( X ,BE ( I A ) , XM ( K  )  ,  0)  ) 

Y=Y+EYY 

13  CONTINUE 

14  I F ( NORM )  17,17,15 

15  DO  16  KA=  1 , NORM 

25  EYY=EA ( I A )*RAMZ(KA ) *01 TON (X,AE{ IA) , RXM ( K A ) ,0 ) - E8 ( I  A ) *RAMZ ( KA  )* 

1  DIT0N(X,BE(IA),RXM(KA),0)  +  ( (EB(IA)-EA(IA) )/{BE( I A ) — A  E  C IA)  ) )  * 
2RAMZI KA)* (DITIN (X ,AE( I  A) ,RXM( KA),0)-  OIT 1 N( X , BE ( IA ) , RXM { KA ) , C )  ) 
Y=Y+EYY 

16  CONTINUE 

17  CONTINUE 
IF(N0K)595,595, 500 

500  DO  59C  I B= 1 ,N0K 

IF  (X  -AG( IB) )590, 590, 509 
509  IF  (NCL) 530, 530,5  10 
SIC  DO  52  C  .1=1,  NCL 

26  EYY=GA( IB)*(QA{  I) *  ENTON ( X, AG ( IB) , A{ I ) , 1 )  +  (  ( QB  ( I  )  -CA(I)  )  / 

1  (2.* (B ( I ) -A ( I  )  ))  )*( ENTON (X, AG{ IB) ,A( I ) ,2 )-ENTON{X , AG ( IB ) , B ( I ) , 2 ) ) 

2-QB ( I ) *ENTON( X, AG { IB) ,B( I ), 1 ) )  +  ( (GB (IB)-GA(IB) )/ CBG( I8)-AG(I8) ) ) 
3*(QA( I)*ENT1N(X,AG( IB), All ),1)  +  ( ( QB ( I ) -CA ( I ) ) / {2 . * ( B ( I ) -A (  I )  )  )  ) 
4*( ENT  IN ( X , AG ( IB ) , A ( I ) , 2 ) -ENT  1 N ( X, AG ( I B ) ,B( I ) ,2 ) )  -  GB ( I ) *ENT IN 
5  (  X ,  AG  (  I B  )  ,  B  (  I  )  ,2  )  ) 

Y=Y-EYY 

27  EYY=-GB( I B ) * ( QA ( I )*ENTCN(X,BG(IB) ,A{ I ), 1 )  + ( ( QB ( I ) -QA ( I ) ) / 

1  (2.*( B( I ) -A ( I ) ) ))*(ENTON(X,BG(IB),A(I ) , 2 )-ENT  ON ( X, BG ( I B) , B ( I ) , 2  )  ) 
2-QB (I  )* ENTON (X,BG ( IB ) ,B( I ) , 1 ) )-  ( (GB( IB)-GA{ IB) )/( 8G(IB)-AG(  IB))) 
3*(QA( I)*ENT1N(X,BG( IB ), All), 1 )+  ( (QB( I)-QA(I) )/(2.*(B(I)-A(I )) )) 
4*(ENT1N(X,BG( IB),  A C I ) ,2)-ENT1N(X, BG( IB) ,B( I) ,2) )  -  QB ( I ) *  ENT  1  N 
5  ( X.BG(IB) ,B( I) ,2  )  ) 

Y=Y-EYY  1 
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520  CONTINUE 

530  IF(N0F)560»560»54C 

540  DO  550  N=1,N0F 

28  EYY=G A  (  I B ) *FY ( N) * ENTON ( X,  AG (  IB)  ,XFO(N)  ,0)  -GB  (  IB  )*FY( N ) »£NTCN ( X, 
1BG( IB)»XFO(N)  ,0)  +  ( (GB( IB)-GA(IB) )  /  < BG ( I B ) -AG ( I  8 )>  )  *  FY ( N  J  * 
2(ENT1N(X,AG( IB)tXFO(N) ,0)  -  ENT  1 N {  X,BG  (  I  B  )  ,  XFO  (  N )  ,  C )  ) 

Y=  Y-E  YY 
550  CONTINUE 
560  IF(NOP) 590,590,570 
570  DO  580  M=1,N0R 

29  EYY  =GA( IB)*RY(M) *ENTCN( X,AG( IB) ,XR(M)  ,0)  -GB { 1 8 ) * RY ( M ) * ENTONI X, 
1 BG( 18 ), XR ( M ) ,0 )  +  ( (GB( IB)-GA(IB) )  /  ( BG ( I B ) - AG ( IB ) ) )  *  RY(M)  * 

2  (  ENT  1N(X,AG(  IB).XR(M)  ,0)  -  ENT  1  N(  X,  BG  (  I B  )  ,  XR  ( f' )  ,0)  ) 

Y=Y-E  Y Y 
580  CONTINUE 
590  CONTINUE 
595  CONTINUE 
RETURN 
6N0 
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FUNCTIGN  FTON(X,A,B,N) 
FTON  =  0. 

I  F  (  X- A  )  3 »  3» 1 

1  IF(X-E)3»3,2 

2  FTON= (X-B )**N 

3  RETURN 
END 


FUNCTION  FT  IN  (X,A,B,N) 

FT  1 N  =0. 

IF  ( X-A ) 3 »  3V 1 

1  IF  ( X-B ) 3 , 3 1 2 

2  FT  1 N  =  ( X— B )** ( N+ 1 )  +  (B-A)*( (X-B)**N) 

3  RETURN 
END 
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FUNCTION  ENTCN  (X,AVB,N) 

FN=N 

ENTON  =  0. 

IF(X-A)5,5,  1 

1  IF(X-e>5,5,2 

2  I F ( A-8  )  3 »  3 » 4 

3  ENTON  =  ( X-B  )  **( N+1)/(FN+1.) 

GO  TO  5 

4  ENTON  =  { X-B )**{N  +  1 ) / { F  N+ 1 . )  -  (  A—  8  )  *  *  (  N  +  1  )  /  (  F  N  ♦  1  . ) 

5  RETURN 
ENO 


FUNCTION  ENTlNt X, A, B,N) 

FN=N 

ENT1N  =  0. 

IF(X-A)5,5, 1 

1  IF(X-B)5»5»2 

2  I F ( A-B ) 3 , 3t 4 

3  ENT1N  »  { X-B )  **( N+2)/(FN  +  2. )  +  ( B-A )*( X-B )** (N+1) /( FN+1 . > 

GO  TO,  5 

4  ENT1N  =  (  X-B )  ■**(  N+2  }  / ( FN+2 .  )  +  (  B-A)  *(X-8)**{N+D  /  CFN+1.) 

I  -(A-B)**(N+2)/(FN+2. >  -(< B-A )*( A-B )**( N+ 1 ))/( FN  +  1  . ) 

5  RETURN 
ENO 
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FUNCTION  DITON(X,A,B,  N) 

FN  =  N 

DITON  =0. 

IF(X-A)5,5f 1 

1  IF(X-E)5,5,2 

2  IF(A-E)3t3»4 

3  DITON  =  (X-B)*»(N+2) /( (FN  +  1 . )*( FN+2. ) ) 

GO  TO  5 

4  DITON  =  (X-B)**(N+2)/((FN+l.)*(FN+2.) )  -  (X*  (  A-9  )  **  (  N+ 1  )  )  /  (  FN  +  1  .  ) 

l-(A-B)**(N+2)/((FN+l.)*(FN+2.))  +  ( A* ( A-B) ** ( N+l )  ) / ( FN+ 1  .  ) 

5  RETURN 
END 


FUNCTION  DIT1N(X ,A,8,N) 

FN=N 

DIT1N=0. 

I F( X-A) 5t  5» 1 

1  I F  (  X-  E  )  5 » 5 1  2 

2  IF(A-E)3,3t4 

3  DITIN  =tX-B)**(N+3) /( (FN+2. )*( FN+3. ) )  +  ( (B-A)*( X-B) *»(N+2 ) 1 / 

1 ( ( FN+  1 . ) *  C  FN+2. ) ) 

GO  TO  5 

4  DITIN  =  (X-B)**(N  +  3)/( (FN  +  2. )*( FN+3. ) )  +  ( ( B-A ) » (X-B 1 ** (N  +  2 )  )  / 

1  (  (FN+1.  )*(  FN+2.  ))-(X*-(  A-B)  *«(N+2>  )  /  (  FN+2  .  )- (  X*  (  B-A  J  *  (  A-B  )  **  (  N+ 1  )  )  / 

2  (FN+1.)  -  (A-B)**(N  +  3)/( (FN+2. )*(FN+3. )  )  -  ( B-A )* ( A-B  )** (N+2 )  / 

3  ( (FN  +  1 . ) * ( FN+2 . )  )  +  (A*(A-B)**(N+2))/(FN+2.)  +  ( A  * ( 8-A ) * ( A-B ) ** 

4  (N+l ))/( FN+1 . ) 

5  RETURN 
END 


I 
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SUBROUTINE  G AUSS2 ( N, E P, A , X, KER 1 
DIMENSION  A ( 22, 23  1  ,X(22) 

NPM=N+1 

10  DO  34  L=  1 » N 
KP*0 

Z  =  0.  C 

00  12  K  =  L  » N 

IF(Z-ABSF(A(K,L)  1)11,12,12 

11  Z-ABSF ( A { K, L ) ) 

K  P=K 

12  CONTINUE 
IF(L-KP) 13,20,20 

13  DO  1H  J=L,NPM 
Z=A(l, J) 

A ( L, J ) ® A ( KP , J 1 

14  A ( KP  ,  J ) =Z 

20  IF(AESFIA(L,L) 1-EP 150,50, 30 

30  I F { L-N 131,40,40  ' 

31  LP 1 =L+ 1 

DO  34  K  =  LP 1 , N 
IF(A(K,L) 132,34, 32 

32  RATIO=A(K,L)/A(L»L) 

DO  32  J=  LP  1 » NPM 

33  A(K, J)=A(K, J)-RATIO*A(L, J) 

34  CONTINUE 

40  DO  42  1=1, N 
1 1 =N+ 1-1 

DO  43  J=l,1 
JPN= J+N 
S=0.0 

I F ( I I-N141 ,43,43 

41  IIP1=II+1 

DO  42  K=IIP1,N 

42  S=S+A( II ,K J*X(K) 

43  X ( 1 1 )  =(A(II,JPNJ-S)/A(II,II) 

KERS  1 

75  RETURN 


6? 


50  K E R=2 
RETURN  . 
ENO 
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